Parallax-Tolerant Image Stitching Based on Robust Elastic Warping

1. Input Images and Basic Setup

2. Detect and Match Features (SIFT)

numPoints = 5000;

 

% SIFT detection. Tune `EdgeThreshold` / `ContrastThreshold` if few points

% appear or too many noisy points appear.

points1 = detectSIFTFeatures(gray1, EdgeThreshold=500, ContrastThreshold=0);

points2 = detectSIFTFeatures(gray2, EdgeThreshold=500, ContrastThreshold=0);

 

% Uniformly sample points (helper function exists in MATLAB File Exchange

% or implemented in this repo). This reduces spatial bias.

points1 = selectUniform(points1, numPoints, size(gray1));

points2 = selectUniform(points2, numPoints, size(gray2));

 

% Extract descriptors and match them

[features1, validPoints1] = extractFeatures(gray1, points1);

[features2, validPoints2] = extractFeatures(gray2, points2);

indexPairs = matchFeatures(features1, features2);

 

% Collect matched point coordinates

matchedPoints1 = validPoints1(indexPairs(:, 1), :);

matchedPoints2 = validPoints2(indexPairs(:, 2), :);

matchedPoints1 = matchedPoints1.Location;

matchedPoints2 = matchedPoints2.Location;

 

% Visual check of matched points (montage view)

figure('Name','SIFT matches');

showMatchedFeatures(img1, img2, matchedPoints1, matchedPoints2, "montage");

title("Points matched using SIFT features")

3. Estimate Global Projective Transform (homography)

4. Create Panorama Reference and Warp By Global Homography

tformInv = invert(tformH);

xLimitsIn = [1, size(img2, 2)];

yLimitsIn = [1, size(img2, 1)];

[xLimitsOutInImg1, yLimitsOutInImg1] = outputLimits(tformInv, xLimitsIn, yLimitsIn);

 

xWorldLimits = [min(xLimitsOutInImg1(1), 1), max(xLimitsOutInImg1(2), size(img1, 2))];

yWorldLimits = [min(yLimitsOutInImg1(1), 1), max(yLimitsOutInImg1(2), size(img1, 1))];

panoWidth = ceil(diff(xWorldLimits));

panoHeight = ceil(diff(yWorldLimits));

panoRef = imref2d([panoHeight, panoWidth], xWorldLimits, yWorldLimits);

 

% Warp both images into the panorama reference (global homography alignment)

warpImg1 = imwarp(img1, rigidtform2d(), outputView=panoRef);

warpImgH2 = imwarp(img2, tformInv, outputView=panoRef);

mask1 = imwarp(ones(size(img1, [1, 2]),'logical'), rigidtform2d(), outputView=panoRef);

mask2 = imwarp(ones(size(img2, [1, 2]),'logical'), tformInv, outputView=panoRef);

 

figure('Name','Global Homography Alignment');

imshowpair(warpImg1, panoRef, warpImgH2, panoRef);

title("Alignment Using Global Homography")

5. Select Control Points and Fit Robust Elastic Warp (REW)

[~, idx1] = unique(round(matchedPoints_ok1), 'rows', 'stable');

[~, idx2] = unique(round(matchedPoints_ok2), 'rows', 'stable');

idx = intersect(idx1, idx2);

 

fixedPtsInImg1 = matchedPoints_ok1(idx, :); % fixed points (Img1 coords)

movingPtsInImg2 = matchedPoints_ok2(idx, :); % corresponding points in Img2

fixedPtsInImg2 = transformPointsForward(tformH, fixedPtsInImg1); % Img2 world coords

 

% Regularization parameter for REW (empirical). You may tune to control

% smoothness vs. data fidelity. Here scaled by image area.

lambda = 0.001 * prod(size(img1, [1, 2]));

 

% Fit REW: this constructs the warp object (obj) and selects inliers.

% The class `RewWarp` is part of this repo and implements the method from

% the paper. Inspect `RewWarp.m` for details.

obj = RewWarp(movingPtsInImg2, fixedPtsInImg2, lambda)

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 2.569684e-18.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.874648e-18.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.015159e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.063735e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.159118e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.227761e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.244940e-17.

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 3.261929e-17.

obj =

RewWarp with properties: movingPts: [165×2 double] fixedPts: [165×2 double] wx: [165×1 double] wy: [165×1 double] ax: 0.0023199 bx: 0.094602 cx: -286.65 ay: 0.0027952 by: 0.003731 cy: -16.556 lambda: 75398 gxy: [165×2 double] inlierIdx: [165×1 double]

% update matched points using 3 sigma removal

fixedPtsInImg1 = fixedPtsInImg1(obj.inlierIdx, :);

movingPtsInImg2 = movingPtsInImg2(obj.inlierIdx, :);

 

figure('Name','Inliers After REW Fit');

showMatchedFeatures(img1, img2, fixedPtsInImg1, movingPtsInImg2, "montage");

title("Inlier matches after deduplication and 3-sigma removal")

6. Prepare Meshgrid Coordinates For Local Warping

u_im = xLimitsOutInImg1(1):xLimitsOutInImg1(2);

v_im = yLimitsOutInImg1(1):yLimitsOutInImg1(2);

[U, V] = meshgrid(u_im, v_im);

[u_im_, v_im_] = transformPointsForward(tformH, U(:), V(:));

 

offset_x = round(xLimitsOutInImg1(1) - panoRef.XWorldLimits(1));

offset_y = round(yLimitsOutInImg1(1) - panoRef.YWorldLimits(1));

 

% Build full panorama grid in world coordinates and transform to image

[u, v] = meshgrid(xWorldLimits(1):xWorldLimits(2), yWorldLimits(1):yWorldLimits(2));

u = imresize(u, [panoHeight, panoWidth]);

v = imresize(v, [panoHeight, panoWidth]);

[u_, v_] = transformPointsForward(tformH, u(:), v(:));

 

% Overlap region between images (used to limit where local warp applies within `warpImage` object-function)

[overlapLimitsU, overlapLimitsV] = outputLimits(tformH, [1, size(img1, 2)], [1, size(img1, 1)]);

 

meshImg2_U = reshape(u_im_, size(U));

meshImg2_V = reshape(v_im_, size(V));

meshPano_U2 = reshape(u_, size(u));

meshPano_V2 = reshape(v_, size(v));

 

% Execute the local REW warp (warpImage is provided with the REW implementation)

[warpImg2, gx, hy, eta] = warpImage(obj, img2, meshImg2_U, meshImg2_V, meshPano_U2, meshPano_V2, offset_x, offset_y, overlapLimitsU, overlapLimitsV);

 

figure('Name','REW Warp Result');

imshowpair(warpImg1, panoRef,warpImg2,panoRef);

title("Alignment Using REW")

7. Blending To Create Final Panorama

maskWarped = warpImg2 > 0;

maskWarped = maskWarped(:, :, 1);

maskOverlap = mask1 & maskWarped;

pano = (im2double(warpImg1) + im2double(warpImg2)) ./ (im2double(maskOverlap) + 1);

 

figure('Name','Final Panorama');

imshow(pano, panoRef);

title("Stitching Using REW")

8. Visualization: heatmap and deformation grid

CAMshow(im2uint8(pano), eta)

title("\eta heatmap")

 

% Visualize deformation grid on original and panorama coordinates

interval_Grid = 40; % grid sampling (pixels) used for visualization

[gridXInImg2, gridYInImg2] = meshgrid(1:interval_Grid:size(img2, 2), 1:interval_Grid:size(img2, 1));

[gridH, gridW] = size(gridXInImg2);

[gridXInWorld, gridYInWorld] = transformPointsInverse(tformH, gridXInImg2(:), gridYInImg2(:));

 

[I, J] = worldToSubscript(panoRef, gridXInWorld, gridYInWorld);

deformXInImg2 = gridXInImg2(:) + diag(gx(I, J));

deformYInImg2 = gridYInImg2(:) + diag(hy(I, J));

[deformGridXInWorld, deformGridYInWorld] = transformPointsInverse(tformH, deformXInImg2(:), deformYInImg2(:));

 

% Show grid on the global homography fused image

globalHomoFuse = (im2double(warpImg1) + im2double(warpImgH2)) ./ (im2double(mask1 & mask2) + 1);

Gridshow(globalHomoFuse, gridXInWorld, gridYInWorld, gridH, gridW, panoRef);

title("Grid Show in Global Homography")

 

% Show grid after adding local TPS-like deformation from REW

Gridshow(pano, deformGridXInWorld, deformGridYInWorld, gridH, gridW, panoRef);

title("Grid Show: Global Homography + Local TPS Warp")

% Show deformation field overlayed on img2

Gridshow(img2, gridXInImg2, gridYInImg2, gridH, gridW, imref2d(size(img2)));

hold on;

quiver(gridXInImg2(:),gridYInImg2(:),diag(gx(I,J)),diag(hy(I,J)),'-r')

title("Grid and Deformation Field in Img2 (deformation vectors)")

9. Create Panorama Reference and Warp By Global Similarity Transformatioin

% Estimate a similarity transform S such that points in img1 map to img2

tformS = estgeotform2d(fixedPtsInImg1, movingPtsInImg2, "similarity");

 

tformInv = invert(tformS);

xLimitsIn = [1, size(img2, 2)];

yLimitsIn = [1, size(img2, 1)];

[xLimitsOutInImg1, yLimitsOutInImg1] = outputLimits(tformInv, xLimitsIn, yLimitsIn);

 

xWorldLimits = [min(xLimitsOutInImg1(1), 1), max(xLimitsOutInImg1(2), size(img1, 2))];

yWorldLimits = [min(yLimitsOutInImg1(1), 1), max(yLimitsOutInImg1(2), size(img1, 1))];

panoWidth = ceil(diff(xWorldLimits));

panoHeight = ceil(diff(yWorldLimits));

panoRef = imref2d([panoHeight, panoWidth], xWorldLimits, yWorldLimits);

 

% Warp both images into the panorama reference (global similarity alignment)

warpImg1 = imwarp(img1, rigidtform2d(), outputView=panoRef);

warpImg2 = imwarp(img2, tformInv, outputView=panoRef);

 

figure('Name','Global Similarity Alignment');

imshowpair(warpImg1, panoRef, warpImg2, panoRef);

title("Alignment Using Global Similarity")

10. Prepare Meshgrid Coordinates For Local Warping With Two Images

To remove the projection distortion from image 2 to image 1, and combined with a global similarity transform, decompose this into two separate homography transformations for the two images,the transformation matrices are ,

As the paper[1] states, we combine the homography and a similarity transformation using the techniques proposed in ANAP; i.e.,

Now we apply for reference image(image1), ,Note that the built-in function transformPointsInverse here performs the coordinate transformation using the inverse of the homography matrix , and (u_, v_) are the coordinates after applying .

% Compute the transformation for the reference image(image1)

[uInImg1,vInImg1] = transformPointsInverse(tformH,u_,v_);

 

% Overlap region between images (used to limit where local warp applies within `warpImage` object-function)

[overlapLimitsU, overlapLimitsV] = outputLimits(tformS, [1, size(img1, 2)], [1, size(img1, 1)]);

 

meshImg2_U = reshape(u_im_, size(U));

meshImg2_V = reshape(v_im_, size(V));

meshPano_U2 = reshape(u_, size(u));

meshPano_V2 = reshape(v_, size(v));

meshPano_U1 = reshape(uInImg1,size(u));

meshPano_V1 = reshape(vInImg1,size(v));

 

% Execute the local REW warp (`warpImage` is provided with the REW implementation)

[warpImg2, gx, hy, eta] = warpImage(obj, img2, meshImg2_U, meshImg2_V, meshPano_U2, meshPano_V2, offset_x, offset_y, overlapLimitsU, overlapLimitsV);

 

% Deformation image1

warpImg1 = imageInterp(img1,meshPano_U1,meshPano_V1);

 

% visualize

figure('Name','REW Warp Result');

imshowpair(warpImg1,panoRef, warpImg2,panoRef);

title("Alignment Using REW With Similarity")

% final panorama

mask1 = imageInterp(ones(size(img1,[1,2])),meshPano_U1,meshPano_V1);

maskWarped = warpImg2 > 0;

mask2 = maskWarped(:, :, 1);

maskOverlap = mask1 & mask2;

pano = (im2double(warpImg1) + im2double(warpImg2)) ./ (im2double(maskOverlap) + 1);

 

figure('Name','Final Panorama');

imshow(pano, panoRef);

title("Stitching Using REW with Similarity")

References

Support Functions (visual helpers)